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Abstract 

Hydrophobic force, interfacial tension, transverse density profile in confined water system are 
addressed from first principles of statistical mechanics in a lattice model for water. Using molec- 
ular mean field theory technique we deduce explicit expressions for each of the above mentioned 
phenomena and show that hydrophobic force is a manifestation of Casimir-like effect due to hy- 
drogen bond fluctuations in confined water. It is largely influenced by the long range correlations 
of orientational fluctuations. All the computations are parameter-free and we compare favorably 
with results of molecular dynamics simulations and experiments. 
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I. INTRODUCTION 



Hydrophobic force causes attraction between non-polar solutes in liquid water and is 
known to be an important driving force in micellar aggregation, cell membrane formation, 
large protein structures [jj. Although hydrophobic force is widely known to be pervading 
distinct physical, chemical and biological phenomena, its occurence as a solvent induced 
interaction was first suggested only in 1945, by Frank and Evans through calorimetric stud- 
ies on hydrocarbons [2] and later elucidated in biological context by Kauzmann 3|. The 
former study noted that transferring a small hydrophobe into water was accompanied by 
unfavorable free energy change g domurafed by entropy reducf.on due to reorganization 
of vicinal water molecules [2[. Hence, two hydrophobes show tendency to coalesce in order 
to minimize the unfavorable free energy. In the following years the applicability of this 
viewpoint for different sizes of hydrophobes such as alkanes, proteins has been discussed 
[sl . Theoretical and simulational investigations were limited to small solutes and interaction 
between them on scales of few Angstroms [6]. In early 1980s the first direct measurement 
of an attractive force between hydrophobic surfaces has been carried out using surface force 
apparatus Large hydrophobic surfaces made up of chemisorbed mica were employed in a 
crossed cylinder geometry and the measured force between them was related to interaction 



a 



The force was seen 



free energy per unit surface area using Derjaguin approximation 
to be influential upto hundreds of Angstroms and stronger than inter-surface van der Waal 
forces j?!. The qualitative aspect of the result i.e., long range and monotonic nature of the 



force withstood the test of time M 



11| . To understand the same, various physical mecha- 



nisms were proposed, some generic to any fluid such as metastability of confined fluid [12| , 



dewetting- induced cavitation under liquid- vapor coexistence conditions [l3|, fluid structur- 
ing effects 



14j | and some, dependent on surface details like correlated dipolar fluctuations 



15|, charged bilayer patches [l6|, nanobubbles [l7|. The above phenomenological descrip- 
tions, however, are either envisaged in narrow range of fluid conditions or specifically depend 
on hydrophobization patterns on surfaces. Besides, they were unsuccessful in reproducing 
generic features of hydrophobic force seen in experiments 9|. We attempt to understand in 
a model study the nature of hydrophobic force by analyzing correlations of hydrogen bond 
fluctuations in liquid water and the effect of modification caused by presence of surfaces on 
these fluctuations. The analysis is carried out at generic fluid conditions within the model 



and applies to generic surfaces. 

We discuss some essential aspects in modelling water and past attempts in this direction. 

n 

Hydrogen bonding is an orient at ion- dependent attraction between two water molecules [18[ . 
In order to analyze hydrogen bond fluctuations it is therefore essential to envisage both 
density and orientational degrees of freedom of each water molecule. There are models 
galore proposed and analyzed to reproduce anomalous thermodynamic properties of wa- 
ter 19|, l20|. However, theoretical attempts to envisage fluctuations in water are limited 
to Ornstein-Zernike-like phenomenological approaches, wherein integral equations only in 
terms of density correlation were heuristically proposed and are numerically solved using 



different closure approximations 2l|, |22| . Wertheim's theory of associating fluids envisages 



similar density correlations to be solved in compliance with steric constraints imposed by 



formation of associated n-mers (23| 
approximation schemes employed [21 



he success of these descriptions crucially depends on 



24l |. Other approaches specific to molecular fluids such 



as RISM were seen to be less predictive in case of associating fluids [25|. Also, perturbation 
theories based on smaflness of interaction strength found little success in reproducing liquid 



phase properties of water-like fluids 



26|. 



Quantum mechanical calculations on water dimer in gas phase and diffraction study of ice 
forms provide sufficient evidence for specific nature of hydrogen bonding j2^, i.e., positively 
polarized hydrogen atom covalently bonded to an oxygen interacts only with negatively po- 
larized lone-pair of neighboring oxygen. The specificity necessitates the density of hydrogen 
bonds and dangling bonds (hydrogens or lone-pairs which are not hydrogen bonded) to be 
commensurate with molecular density, stated as a sum rule 28|. Consequently, density and 
orientational fiuctuations (the latter being inherently connected to bond fiuctuations) are 
not totally independent; their long wavelength fiuctuations, especially, are to be consistent 
with the sum rule. Effective interaction models for water designed for numerical simula- 
tions (molecular dynamics) 20|, l29| provide successful instances of (implictly) envisaging 



most essential features of hydrogen bond fiuctuations consistent with the sum rule. A water 
molecule is often modelled as a polar molecule with charges corresponding to hydrogens and 
lone-pairs placed at vertices of a tetrahedron 29|, |30| • A complete description of molecular 
correlations in such models can be achieved by defining a set of orthonormal vectors in terms 
of atomic coordinates and defining correlations among them. Large scale molecular dynam- 
ics (MD) simulations at ambient conditions reveal that density correlations are short ranged 



and extend no further than 12 A; whereas, dipolar orientations of water molecule, which 
are receptive to bond fluctuations in the neighborhood, are correlated over large distances, 
atleast upto 75 A [31]. Two correlation lengths of order 5.2 A and 24 A were inferred. 
Coulomb interactions, surprisingly, have little effect on asymptotic behavior of these cor- 



relations [311 • It is also suggested that mesoscopic hydrophobic solutes when coupled to 



dipolar fluctuations in water medium experience a long range exponential force, which is 



attractive in nature and dependent on shape and mutual orientation of solute surfaces |31|. 
The exponential decay bears a striking consistency with that seen in experiments measuring 
hydrophobic force between large surfaces j3, 0] . For the case of large hydrophobic surfaces, 
correlations in confined water need to be ascertained. To simulate such a system the surfaces 
need to be several times larger than the longest correlation length in the system in order 
to obtain proper equilibrium conditions. This requires huge system size that would render 
the simulation prohibitively resource intensive. In addition, the accompanying free energy 
change could be very small due to weak nature of correlations at large distances. Instead, 
we take analytic route to describe hydrogen bond fluctuations in water and the effect of 
spatial confinement on them. 

Hydrophobes are known to interact unfavorably with water molecules in contrast to strong 
and cohesive water-water interactions. In addition, large surfaces substantially disrupt the 
hydrogen bond network whose fluctuations are suppressed at surface boundaries. The setting 
is ideally suited for fluctuation-induced force between the surfaces driven by thermal energy 
in confined water. Forces of this nature are generically called Casimir forces as they were first 
discussed by Casimir in the case of electromagnetic fluctuations confined between conducting 



plates [32| which was later studied in detail by Lifshitz j33|]. Such forces are now envisaged 
in widely different contexts 3^]. Fisher and de Gennes argued that when a binary liquid 
mixture is confined between surfaces which have specific affinity towards one of the fluid 
components, Casimir-like density fluctuations in the liquid give rise to an effective force 



between surfaces 35|. Origin of the force is entropic in nature; in that, the free energy 
is increased due to restriction imposed on fluctuations by the boundaries, thereby system 
tends to minimize the separation in order to reduce the free energy cost. We investigate 
hydrophobic force to be a manifestation of Casimir-like force, here, due to density and 
orientational fluctuations in liquid water. 

When water is conflned between hydrophobic surfaces the inherent fleld fluctuations van- 



ish on surfaces. Furthermore, water molecules at interface with each surface have restricted 
orientational entropy owing to repulsive hydrophobe-water interactions. This effect gives 
rise to modified fluctuations at each interface. We study the collective consequences of these 
effects on the overall free energy of the system in a model study. We define a simple water 
model wherein density and orientations of a water molecule are envisaged. The specific 
nature of hydrogen bonding between molecules is incorporated and associated constraints 
on the bond network are taken care of in the analysis. Partition function is solved about a 
mean field which is consistently deduced at arbitrary densities within the model. Fluctuation 
properties are also deduced. Large correlation length is seen for orientational fluctuations. 
Two macroscopic surfaces are envisaged as boundaries in a spatial dimension. The change 
in free energy in the system due to the presence of surfaces is calculated and is seen to 
be composed of three important contributions : (i) Casimir part, which arises solely from 
discretization of fluctuation modes between boundaries and is generic to all surfaces; (ii) 
Interfacial tension, which is free energy change due to modified fluctuations at hydrophobe- 
water interface. It is dependent on nature of surface-water interaction and to a small extent, 
also on separation distance between the surfaces; (iii) Interfacial fluctuations-induced free 
energy, which is due to correlations of modified fluctuations at both interfaces. It depends 
on type of both surfaces and their interaction with water. The results are discussed for dif- 
ferent types of surfaces such as hydrophobic and hydrophilic. We flnd that the Casimir part 
is leading contribution and is an inverse power-law function of separation distance. However, 
numerically the magnitude of Casimir part is signiflcant only upto four times the longest 
correlation length in the model. Interfacial tension also varies with separation distance, 
but its variation is numerically insigniflcant. Interfacial fluctuations- induced contribution is 
seen to be exponential 



y decaying with distance, analogous to the force form deduced for 
mesoscopic surfaces 3l|- Furthermore, we flnd that all the contributions are of comparable 
order of magnitude consistent with experimental values. The dependence of the force on 
fluid conditions like temperature, average hydrogen bonds is also discussed. Our results 
indicate that hydrophobic force qualitatively imitates Casimir-like force behavior within our 
model study. It is desirable to emulate the computation within more realistic models of 
water possibly with the help of MD simulations. We also looked at transverse density proflle 
for conflned water and show that an increase in density occurs near interfaces. 
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II. MODEL FOR WATER 



We envisage our water model in the background of a lattice to exactly account for hard- 
sphere repulsion between atoms at short distances. Also, constraints of hydrogen bond 



network are explicitly taken care of in lattic^ 
field (MMF) technique in an earlier paper 



28 



Dackground. We developed molecular mean 
to address the same model in the infinite 



volume limit. We present here a brief summary of the model and MMF technique and then, 
address the case of confined water. 

On a three dimensional hypercubic lattice we define occupation field W{r) = {0, 1} 
corresponding to water being absent or present, respectively, at the site r = {x,y,z). At 
each occupied site we define bond arm field Ha{r) which resides on the links around the site 
r. Ha{r) = {0, ±1}, where a = {±1,±2,±3} correspond to six directions around the site. 
Ha{i^) = 1 refers to hydrogen arm on the corresponding link, —1 to lone pair arm and for 
no arm. The constraints between W{r) and Ha{r) being, 

J^Hlir) = mr) (1) 

a 

J2Ha{r) = (2) 

a 

which imply that every water molecule has two hydrogen and two lone-pair arms only. A 
hydrogen bond is realized when two water molecules two lattice units apart have one of 
each's hydrogen and lone-pair arms meet at a site, as shown in Fig.([T]). When two molecules 
are on near- neighbor sites they are disallowed to have any non-zero bond arm on the link 
between them. The constraint is given by : 

W{r) (^H^ir + eo.)^ =0 (3) 

We write a general interaction Hamiltonian in terms of field as below : 

' a, a 

where, A is interaction strength and a, a denote directions around site r. There are addi- 
tional restrictions on field, namely, 
(i) at any site no more than two bond arms meet i.e., 

0<yHlir + e^)<2 (5) 



a 



(ii) two non-zero bond arms are disallowed from meeting at a site i.e., anti-bonds are disal- 
lowed, 

-1 < 5^i^a(r + e,) < 1 (6) 

a 

The grand canonical partition function for the system at a finite chemical potential jl for 
water and inverse temperature (3 is given by : 



n n 



(7) 



where, the prime over summation indicates that the W{r) and Ha{r) sum have to be carried 
out in compliance with Eqs.( TTP|5|5p]) . Evaluating Z amounts to enumerating all possible 
bond configurations that satisfy above constraints and calculating the exponential in Eq.(I7]) 
for those configurations over allowed range of W and Ha at fixed values of /i, (3 and V the 
volume of the system. 

The restrictions represented by Eqs. fl3]|5|6|) are at sites where there is no water. These 
are shown in Fig.(l2]). To implement them in our analysis it is useful to define two discrete 
integer fields 6(r), g(r) : 

b{r) = Y,Hl{r + ea) (8a) 

a 

q{r) = Y,Ha{r + eo) (8b) 

a 

The discrete field 6(r) counts the number of non-zero arms in the neighborhood of site r, 
while q{r) measures the charge i.e, difference between number of hydrogen arms and lone-pair 
arms meeting at site r. By construction, 6(r) varies between and 6 on a three dimensional 
hyper-cubic lattice and g(r) in turn varies between —h{r) to 6(r). By imposing the condition 
that h{r) < 2 in our analysis we ensured that no more than two arms can meet at a site. 
Furthermore, for 6(r) = 2 we demand g(r) = to disallow anti-bond configurations. In 
terms of these variables, Eqs.( l3|5|6|) can be rewritten as : 

W{r)h{r) = (9) 
(6(r),g(r)) = {(0, 0), (1, 1), (1, -1), (2, 0)} (10) 

where, b{r) and g(r) values are restricted only to the above set of mutually exclusive pairs. 
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We now rewrite the partition function as below : 



n E 

^ W{r),H^{r) 
b{r),q{r) 



-/3^(H-z>g2(r)-/iiy(r)) 



(11) 



where, we have additionally introduced a chemical potential u for dangling bond configura- 
tion i.e., (6, q) = (1, ±1). The fields b, q are summed over allowed range given in Eq. fllU]) and 
the prime over summation indicates that Eqs.( TT|2|5|5]) act as constraints in the evaluation. 
Note that, since only hydrogen bond interaction is envisaged in the model, the Hamiltonian 
7i can be rewritten as a simple expression : 



(12) 



where, Kronecker delta function denoted here as 6{p, q) is defined as 5{p, q) = 1 for p = q 
and otherwise. All the possible hydrogen bond configurations are implied from solving 
the non-local constraints Eq.(|S]). These constraints are enforced in the partition function by 
introducing auxiliary fields, as given below : 

\ a / ^(^^-^^ 



exp 



TT 

'N 



Vir) lbir)-J2Hl{r + e^] 



6 lq{r),J2Ha{r + ea] 



2N 



(13a) 
(13b) 



where, //(r) and 0(r) act as dual fields to density and charge of bond arms in a local 
neighborhood. The discrete i] and fields take integer values in the range [— A^, A^] at every 
site, where A^ is any suitably large integer (greater than 8). 

The partition function can be rewritten in terms of new variables and auxiliary fields as : 



n E 



E 



(2Ar + l)2 
W{r),Ha{r) ^(r),0(r) 
b{r),q{r) 



exp 



E 



P{n - uq'{r) - jlW{r)) 



N' 



(14) 



Here, prime over summation refers to sum being restricted to local on-site constraints 
Eqs.( IT][2][9|) only. The introduction of auxiliary fields ri{r) and 0(r) allows summation over 
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other discrete fields W{r), Hc,{r), b{r), q{r) within their respective allowed range at each site 
without any restriction from the neighborhood configurations i.e., as if a single site functional 
Zsite, as given below : 

Z= /"MMn^..te(r/(r),0(r),V„r/,V,0) (15) 

r 

where, the summation over rj, (f) fields is transformed into an integral in the limit of — t- cxo 
and Zaite is given as sum of weights corresponding to each allowed state i.e., void, dangling 
bond, hydrogen bond and water. It is given by, 

Zaue = 1 + 2z/e-*''('') cos(0(r)) + Ae'^^''^'^) + ^^C{r^, 0, V^r^, V„0)) (16) 

C(?7, V„?7, Va</)) = ^ exp 

= 0, ±1 
a = ±1,±2,±3 

where, u = exp(/3z>), A = exp(/3A) and /i = exp(/3/i) are fugacities of dangling bond, hydrogen 
bond and water states, respectively. The orientational degrees of freedom of water yields 
C{r], 0, VaV^ Vq,0) given by Eg. (1171) . where the summation is over orientations at site r. The 
prime over summation indicates H^s of each orientation satisfy Eqs. (fT][2!) . The exponential 
corresponds to an orientation and it is a function of dual fields at near-neighbor sites towards 
which non-zero bond arms of the orientation are directed. The densities of dangling bond 
(DB), hydrogen bond (HB) and water (p) are calculated from partial derivative of partition 
function with respect to /3z/, f3\, /3fi, respectively. 



J2iHl{r)7]{r + e„) + H^{r)(t){r + e„)) 



:i7) 



A. MMF theory 

The partition function has a unique maximum at isotropic and homogeneous field con- 
figuration = = 0. Zsite at the maximum is given by Zo : 

Zo = (l + 2i/ + A + 90/x) (18) 

The extremization condition 'Ylixi^ I ^^)^ ~ ^ implies the sum rule of the system i.e., DB + 



2HB = 4p exactly 28|], while extremization with respect to field is trivially satisfied. To 
the zeroth order, partition function is Z = {Zo)^ and sum rule translates as : 

2i^ + 2A = 4(90p) (19) 



Using this relation, the densities of danghng bond, hydrogen bond and water are given upto 
zeroth order as : 

DB = iy-^(\nZ) = — (20a) 

Ehminating A from equations for DB, HB we obtain : 

HB = 2p--^(l-3p) (21) 

We call Eq. fl2T]) the equation of network. It is a manifestation of sum rule in terms of model 
parameters. We choose dangling bond energy parameter to be zero i.e. u = 1 and measure 
temperature (/3~^) in units of hydrogen bond strength (A). To zeroth order the theory is 
now parameter-free and all densities can be obtained as a function of temperature only. The 
equation of network can also be written in terms of average hydrogen bonds per molecule 

h=^ as: 

p 

h = 7-- (22) 

P 

From zeroth order partition function the mean field free energy Gm per unit volume can be 
given in terms of densities as : 

/3G^ = ln(l - 5p + HB) (23) 

Eq. ( l23i) is analogous to equation of state for a system of hard spheres. It correctly predicts 
the density saturation in the model at p = |, HB = | or /i = 4. Thus, equation of network 
is a manifestation and density saturation effect is a direct consequence of the sum rule. 

Using Eq. f|T6l) for ZgHe we expand dual fields upto quadratic order about their maximum, 
perform Fourier transform on their fluctuations in a large cubic box using periodic boundary 
conditions, then integrate the resulting Gaussian functional in Eq.f llSI) over all field config- 
urations. This yields total free energy per unit volume G upto one-loop correction which 
includes leading contributions due to fluctuations in density and orientations : 

TT 

/3G = + - y ^In {p^^{k)P^^{k)) (24) 
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where, 



P^,(fc) = 64p 



96p 



9 
10 



A(l-A) 



A 



20(|t - p) 



+ 



3(i-p) 



5(DB) 
96p 



(25a) 
(25b) 



and A = 

'^^(1 — cos(fcj)); fcj are vector components of k. Prjn, -P,/,,/, are fluctuation 
propagators of dual fields. The above expressions for propagators are simplified to the 
leading order using sum rule. Precise expressions in terms of original fugacities are given in 
Appendix (jAj). The term involving propagators in the free energy expression (Eq j2^ is the 
entropy contribution about the mean field. 

The correlation functions in the system in the momentum space are given by propagators 
Prjn and P,^,^. Density correlations are dominated by //(r) correlations. They display coordi- 
nation peaks in position space reminiscent of radial distribution function of fluids 2l| and 
do not have any long distance behavior (FigHj). Orientational correlations are dominated 
by 0(r) correlations at large distances. They display a correlation length of upto 4 lattice 
units in liquid phase (FigJl]). 

We also envisaged Coulomb interaction between bond arm charges via a new dual field 
that couples to the charges. We find that Coulomb interactions have little effect on the 
asymptotic behavior of orientational fluctuations and also on MMF results like equation of 



network and equation of state 



28|. 



III. WATER CONFINED BETWEEN MACROSCOPIC SURFACES 

We now study the case of water confined between two macroscopic hydrophobic surfaces. 
As a result of confinement the structure of fluctuations in the system is restricted by the 
boundaries, thereby causing entropy to be a function of separation distance between surfaces. 
In analogy with Casimir interaction, the distance-dependent entropy component of free 
energy of conflned water leads to an effective interaction between hydrophobic surfaces. In 
addition, due to hydrophobe-water interactions, orientational fluctuations are modifled at 
interface of each surface. The modifled fluctuations and their correlations lead to interfacial 
tension proportional to area of surface and an induced interaction between the surfaces. 
The net effect is an interaction force that acts over distances longer than typical hydration 

11 



structure of water. We utilize MMF framework to analyze these effects in a unified fashion 
within the proposed water model. 

We envisage surfaces in the (x, y) plane of rectangular coordinate system; one present at 
z = and other a,t z = L (FigJS])- Each surface excludes water from its region of occupation. 
Hence, W = Q on surface sites. On the immediate layer, i.e., aXz=loiz = L — 1 called 
interface layer, water can be present and can take various orientations. For a hydrophobic 
surface if a non-zero bond arm of interface water is directed towards the surface, there would 
be a dangling bond on surface site; else a void state occurs. There can never be a hydrogen 
bond on surface i.e., 6 7^ 2 on surface. We will take care of these possibilities explicitly in 
our analysis. Consequently, we need not introduce r] and integrals (Eq lT^ on the surface. 
Alternatively, we set = = on surfaces. 

The calculation of partition function begins with formulating the site functional ZsUe 
at each site, which comprises weights corresponding to each allowed state in the model. 
The site functional for all the sites in bulk region is of same form as given by Eq.( !T6|) . 
On interface sites, weights corresponding to void state, dangling bond and hydrogen bond 
states remain unaltered. When a water molecule is present on a interface site its bond arms 
can orient in all possible ways. Only if one of the arms is towards the surface we assign 
a weight exp(/3i/s) to the corresponding orientation. For an ideal hydrophobic surface i.e., 
which is indifferent to bond arms of vicinal water, z/g = (in general, z/5 can be positive or 
negative). Consequently, orientational weights for a water state on interface (with surface 
in 63 direction) are given by : 

/ \ 



interface 



E 



+ exp(/3z/5.) 



E 



a ^ 3 if« = 0, ±1 

iJ3 = 



a ^ 3 i7« = 0,±1 
H^ = ±l 



exp 



; ^{Hl{r)r]{r + e„) + H^{r)<P{r + e„)) 



C(r/, 0) + i^s 



E 



exp 



a ^ 3 /f^ = 0, ±1 



; ^{Hl{r)7]{r + e„) + Ha{r)(j){r + e„)) 



^ CivA) + i^sC'ivA) (26) 
where, prime over Ha sum implies constraints Eqs. (|Tp|) . C'{r],(p) corresponds to affected 



12 



orientations only i.e., those with i^s = ±1 and, z/5 = exp^Sus) — 1 is a function of surface- 
water interaction strength. The site functional Zj for any interfacial site can be arranged as 

Zi = Zsite + i^sfi'C'ivA) (27) 

By definition, us ranges from —1 to 00. We remark that for a perfect hydrophobic surface, 
us = 0. 

The modified site functional at interfacial sites can be recast in the expression for full 
partition function, such that the following decomposition is deduced: 

mmU^sUe n (l + r(ri)) n (l + r(r2))) 

= Z /exp I ^ ln(l + Tin)) + J2 + r(r2)) J \ (28) 

where, Z\\ is partition function for the system with surfaces; Z is for corresponding unper- 
turbed case {us = 0) with 77 = = on surfaces and r(r) is defined only on interfacial sites. 
It is relative orientational weight of affected orientations with respect to Zgitg, i.e.. 

The partition function for unperturbed case Z can be evaluated using MMF technique. 
The leading mean field energy is obtained from the maximum of Zgite at each site and 
fluctuations in rj, fields are analyzed subject to vanishing boundary conditions on the 
surfaces. The interfaces-dependent part in Z\\ is evaluated using cluster technique and the 
corresponding free energy is obtained. The resulting form of total free energy Gtot per unit 
lattice area is organized to be: 

Gtot = Go + Gc + 7si + 752 + Gr (30) 

where. Go + Gc is the free energy obtained from evaluation of Z, analogous to Eq.( l24|) . 
Go includes leading terms proportional to L and constants obtained in large L limit. They 
contribute only to bulk pressure of the system. Gc is the remaining L-dependent part. 7^^, 
7S2 are free energy contributions due to surface-water interaction and evaluated only on 
respective interface sites Ji and I2 respectively. Gr constitutes terms which involve sites of 
both interfaces. Expression for each of the terms is deduced in the remaining section and 
their relevance to hydrophobic interaction is elucidated. 
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We first evaluate Z using the MMF technique described in the previous section. We 
identify the maximum of the functional to be at r/ = = 0. It yields mean field free energy, 
which to the leading order is given by LGm (Eq J2^ . The dual fields are then expanded upto 
quadratic order about their maximum and the resulting Gaussian functional is integrated 
over all field configurations to give one-loop contribution to free energy. In the process, the 
following Fourier transform is employed which satisfies the boundary conditions : 

r]{r) = — — N9 ^ fi{k) exp{ikix + ik2y) sin {mrz/L) (31) 

L J (27r)^ 

— TT 

where r = (x, y, z) is position vector and k = {ki, k2, k^ = ^) denote modes in momentum 
space. Similarly for field. 

The entropy contribution to free energy for unperturbed system is discrete analog of that 
of bulk water (Eq J2^ . in that the integral over wavevector in z-direction is replaced by a 
summation over a restricted number of wavevectors i.e., k^ = ^, . . . , ZllLiil. Xo analyze 
L-dependence, we define entropy contribution per unit area in each mode in z-direction as : 

TT 

Siks) = \ I ^^^^In {P,,Ck)P,m) (32) 

where, the propagators Pr^q, P^^ are the same as in the case of bulk water. Total entropy 
contribution to free energy of confined water is S{k^) summed over allowed values of k^. Its 
large-L behavior can be enumerated using Euler-Maclaurin series expansion [sgI] : 

f(L-l) ^ 

J2 Sik,) = L I ^S{k,) - \ (S(0) + 5(7r)) + (33) 

'=3 = f 

On right hand side of Eq. (l33l) . first term is total entropy contribution in the same volume of 
bulk water. S'(O) and ^'(vr) are free energy densities in modes /cs = and k^ = it respectively. 
They are independent of L. From Eq.f lSSp we infer Gc to be the net difference in entropy 
contribution per unit area between confined water and bulk water in the same volume. Gc 
can be calculated as a series expansion in the leading term being : 

d 



IT 

C 



I3G- ~ 



BoL 



d 

S{k,) 



dk 



-Sik, 



for large L (34) 



fc3=0_ 

where, i?2 = 2 is first Bernoulli constant. Gc is analogous to the Casimir interaction energy 



derived for the case of conducting plates confining electromagnetic fluctuations 



32l |. Hence, 



we call Gc the Casimir part of free energy. It falls-off asymptotically as for large L. 
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In the expression for partition function (Eq j28|) average over interface terms is now pur- 
sued. At each interfacial site, ln(l + r(r)) ~ r(r) is the leading order term. This is justified 
because in Eq. 



for r(r), we note that ^'^ ^^''^'^ 



less than 1, since p < \ and 



90 



90 

< |, C'(0,0) = 60. From Eq. 



whose maximum value is always 



the leading order 



contribution from interface terms is given by : 



§ = (exp($^r(n)+$^r(r2)) 

\ Vne/i r2e/2 / 



(35) 



The average can be evaluated using cluster technique [37j]. Terms that involve sites of same 
interface and those involving sites of both interfaces are segregated. 7s is defined to consti- 
tute terms corresponding to sites on same interface. Each of them is proportional to us or 
its higher order. 75 is given to the leading order as : 




E r(ri)r(r2) 




(36) 



where, A is area of the surface. 75 arises due to surface-water interaction and consequent 
effect on orientational fluctuations in the interfacial region. 

Each of the averages in Eq. fl36l) can be evaluated using a functional integration relation jsS] . 
For an interface site with surface in 63 direction, using Eqs. fl26ll29p (r(r)) is given to the 
leading order as : 



(r(r)) 



H3 = ±l 



J2{Hl{r)Hl,{r)g,ir + eo^,r + e^>) 



(37) 



where, the Ha summation is over affected orientations at site r. The prime over summation 
indicates i^a's of each orientation satisfy Eqs. (fT]l2|) . The exponential in Eq. (l37|) corresponds 
to one such orientation, if^, H^i are bond arms of the same orientation; r + e^, "f" ^ 
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the bond arm locations. The average (r(ri)r(r2)) is given to leading order as : 



(r(n)r(r.)) = ^ ^ 

Hs = ±1 Hs = ±1 

a,a' 



+ {Hl{n)Hl{r2)g^{n + e„, ra + e,) + H^{n)H,{r2)g^{n + e„, rs + e,)) 



(38) 



where, Ha, H^i are bond arms of an affected orientation at site ri and H,^, H^i are those 
of an orientation at site r2. The exponential corresponds to product of the two orientations 
and the summation is over all possible products. The two-point Green's function Qn{Ti-,r2) 
for ?7-field fluctuations between any two arbitrary sites ri = (xi, yi, Zi) and r2 = (x2, 2/2, ^2) 
is given by, 



T 1 " 

2 /■ 

Qr,{ri,r2) = tJ2 

n=l 



{dki){dk2) ,. N •; ^ sin (n7r2;i/L) sin (n7rz2/-^^) 



{2ny 



exp - X2) + ik2{yi - ^2)) 



PrfTjik) 



(39) 



Similarly, ^,/,(ri,r2) for field can be defined using the propagator P^^{k). 

The expression for 75 indicates that it varies with separation distance, owing to the 
L-dependent Green's functions. The asymptotic value of 75 is the interfacial tension for 
hydrophobic surface in contact with water. The leading correction term is proportional to 

for large-L and contributes to force between the surfaces. 

From the cluster expansion of partition function, terms that involve sites of both interfaces 
are grouped as Gr- It is given to the leading order as : 



-(5GtA 



(40) 



Effectively, Gr is connected correlation between orientational fluctuations of both interfaces. 
Hence, we call this contribution interfacial fluctuations-induced part of free energy. The 
averages in Eq. (l40|) can be evaluated using Eq.( l37|l with us corresponding to each interface 
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and using Eq.( l38|) with proportionality factor {1^81^82) instead of {i^sY- The identity of sites 
is as per given in the expression for Gr (Eq j40|) . 



The long distance behavior of Gr is dominated by 0(r) correlations, T]{r) being short 
ranged. Between two hydrophobic surfaces, to the leading order Gr is proportional to 
square of orientational correlations i.e., (^(/.(t"))^, where ^^(r) is an exponentially falling-off 
function for large r (Appendix R|) . 

For the case of mesoscopic surfaces hydrophobic force is suggested to arise from orien- 
tational correlations between water molecules at both interfaces 3l||- The force is seen to 
decay exponentially with separation distance, asymptotically. Gr is thus analogous to hy- 
drophobic interaction free energy of mesoscopic surfaces. However for macroscopic surfaces, 
in addition to Gr, hydrophobic force obtains contributions from Casimir part and interfa- 
cial tension. This distinguishes hydrophobic interaction between large surfaces from that 
of between small surfaces both qualitatively and quantitatively. The non-additive nature of 
lydrophobic interaction with increasing size of surfaces has attracted considerable attention 
3^ and our work provides a direction to elucidate the size dependence in terms of hydrogen 
bond fluctuations in water. 



A. Hydrophilic surfaces 

We can envisage surfaces of generic heterogeniety in our calculation. The heterogeniety 
could be in terms of space-dependent us and/or charge on surface. One of the simplest cases 
is a homogeneous hydrophilic surface with a fixed charge on each site. We first consider 
the case of a positively charged hydrophilic surface. On its interface, the site functional 
comprises weights corresponding to all states. When a water molecule is present on interface, 
its hydrogen arm is restricted from pointing in surface direction. We assign an energetic 
penalty to such orientations and the site functional can be arranged, analogous to the case 
of a hydrophobic surface, as : 



Zj = Z,ite + vsiiC'{ri,(t)) 
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Here, z/5 G (—1, 0) (ideally, vs = —1) and the orientational weights corresponding to affected 
orientations C' {rf^cj)) are given by : 



E 



exp 



a ^ 3 = 0, ±1 



J2iHl{r)r]{r + e«) + H^{r)<j){r + e«)) 



(41) 



The above expression is for an interface site with surface in 63 direction. A negatively 
charged hydrophilic surface can also be envisaged such that for interface water orientations 
with lone-pair arm in surface direction are energetically penalized. Here, the weights for 
affected orientations are : 



E 



exp 



a ^ 3 = 0, ±1 



Y^{Hl{r)r]{r + e„) + H^{r)(f){r + e„)) 



(42) 



We now compute the free energy components Gc, 75i, 752) using their respective 
expressions for different types of surfaces, us is an arbitrary parameter in the calculation. 
It is chosen close to its ideal value for each surface type. The properties of water enter 
the computation via Green's functions Q^f,. These are computed within the model using 
Eq.( l39|) . Due to L-dependent modes in the confined direction, all the free energy components 
that depend on fluctuations are expected to vary with separation distance L. 



IV. RESULTS : HYDROPHOBIC FORCE, INTERFACIAL TENSION 

We first mention that this computation is totally parameter-free on lattice. Hence, the 
best way to interpret results is in terms of physically observable quantities such as p and 
hydrogen bond density. Indeed because of equation of network (Eq j2T]) only one of them is 
independent. We find that it is best to describe in terms oi h = the average number 
of hydrogen bonds per molecule. Temperature is conjugate to HB (total hydrogen bond 
density) and hence it is also implicitly fixed self-consistently due to equation of network, 
as shown in Fig.(jl]). The relation between p and h is simple at zeroth order in MMF 
theory (Eql221), but it becomes non-linear at one-loop level. (Zeroth order is still a reason- 
able app— ion MO Hence. aU de.it.e. DB, HB, , a„d the .ee ene... _n. 
given by Eqs. fl33f36f40l) are evaluated from partition function upto one-loop order using the 
corresponding expressions for propagators (Appendix R|) . 
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In our model MMF theory describes liquid for h > 3 reasonably consistently. For 
2 < h < 3 MMF approximation is not seen to be good, namely, one-loop order terms 
are either comparable or exceed zero-loop term. So, we choose to present our results for 
h > 3. All the potentials and energies are computed in the units of hydrogen bond strength 
A taken to be unity. The lattice constant in the model is arbitrary. By computing physi- 
cal lengthscales such as correlation length it can be fixed. Correlation lengths for density 
and orientational fluctuations (^^) to the leading order are simple expressions given in 
Appendix ([A]), but a precise expression to one- loop order is implicitly given. In Fig.(|l]) we 
plot correlation lengths as a function of h. is only about one lattice unit in liquid phase 
and does not vary considerably with h, while increases with h. In MD simulation density 
correlation length is not seen; this is consistent with MMF result since is equal to the 
minimum length possible in the model and also independent of h. Orientational correlation 
lengths inferred from MD simulation are 5.2 A and 24 A, of which the latter is weaker in 



strength (one-tenth) relative to the shorter one 3l|]- In our water model we have only one 



orientational correlation length which we relate to 5.2 A. For liquid water h value is 



suggested to be about 3.6 (30[. From Fig.(jlj) h = 3.58 corresponds to ~ 3.3 lattice units. 
Consequently, we infer that 1 lattice unit ~ |^ = 1.57 A. 

In Fig. ([5]) various contributions to interaction free energy and their relative magnitudes 
are plotted as a function of separation distance L between surfaces. The plot is presented 
for h = 3.58. Casimir part Gc gives the most attractive force, followed by Gr, while the 
interfacial term 75 is repulsive, albeit very small. Gc, 'Js fall-off as for large L from our 
analytic calculations. Numerically, beyond 15 lattice units they are insignificant. All the 
plots are presented for lattice distance L > 5. For smaller L the results are predominantly 
influenced by surface effects. In the model, for L = 4 there is only one layer which can have 
free orientations (besides two interface layers), while for L > 5 there are two or more such 
layers. 

Force is computed as discrete derivative of total free energy with respect to L and plotted 
in Fig. ([6]) for various h. The curves effectively show that the force can manifest upto a length 
of about 15 lattice units which translates to about four times the orientational correlation 
length in the model. All the free energy components and force obtain major contributions 
from orientational fluctuations. 

Figs.( I7f8|9|l display the /i-dependence of Gc, 75 and Gr functions. The Casimir part 

19 



Gc monotonically increases in magnitude with h. The interfacial tension contribution 75 
decreases and is always slightly repulsive. Interfacial fluctuations-induced part Gy increases 
with h for L > 6. At shorter distances it decreases with increasing h. This indicates that 
the adhesion strength of Gy component is higher for higher temperatures. This behavior 
is qualitatively similar to the temperature dependence of interaction free energy for meso- 
scopic hydrophobic surfaces {2!, H]. This reaffirms our interpretation that Gy component is 
analogous to hydrophobic interaction free energy for mesoscopic surfaces. 

Fig.( !TO|) is the plot for Gr contribution between two hydrophilic surfaces, both of same 
type (hydrogen donor /acceptor) and of dissimilar type. Gr in this case is proportional to 

and hence, correlation length is twice as longer in range than in the case of hydrophobic 
surfaces (where Gr is proportional to (^<^)^). At short distances it is seen to be attractive 
for both combinations. However, for large distances it is weakly repulsive between like- 
charged surfaces, in contrast to attraction between oppositely charged surfaces. Fig.( !TT|) 
depicts the force between hydrophilic surfaces for both similar and dissimilar combinations. 
As expected, the dissimilar pair of surfaces have marginally larger attraction than that of 
similar surfaces. It is interesting to note that like-charged hydrophilic surfaces also have a 
net attraction. This is due to dominance of Casimir part Gc which is indifferent to surface 
charge. 

Fig.f lT2|) displays force between a hydrophobic and hydrophilic surface. It bears similar 
profile as in the case of two hydrophobic surfaces. This is expected because essentially Gy 
is qualitatively same for both cases i.e., proportional to {G,),Y- For all surface combinations 
the force is seen to increase in magnitude with dominantly due to indifference of Casimir 
part to surface types. This is a consequence of the fact that the entropy induced forces are 
largely charge neutral. 

Next, we make an attempt to relate our computational results to those of experiments. 
The free energy values presented in the graphs are in the units where hydrogen bond strength 
is unity. Generally, dimensionful quantities in lattice models and those in corresponding con- 
tinuum models are not the same. So it is best to compare dimensionless quantities. In our 

instance, for h = 3.58 and L = Q lattice units which translates to 6 x 1.57A ^ 9.5A, 

\GtoM - Gtotioo)\ 9x10-^ _„2 ^ .... 

— ~ 10 . rrom experiments mteraction free energy esti- 

7s(cx)) 8.5 X 10-3 I — I 

mate when two hydrophobic plates are about 10 A apart is about 1 mJ m~^ 11 1, while 



interfacial tension is in the range 50 — 100 mJ m 
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4l|; their ratio agrees with our compu- 



tation. In experiments the free energy values are also measured for larger distances all the 
way up to 100 A. Unfortunately our model is not good for these distances. This discrep- 
ancy was already noticed when our results were compared with MD simulation. The simple 
water model has only one orientational correlation length, while there are more than one in 
both MD simulations 28|] and surface force apparatus experiments We conclude that 
while order of magnitude estimate of the strength of hydrophobic force is in agreement with 
Casimir-like energies envisaged here, a few more important details are perhaps missing in 
our simple model of water. 

V. TRANSVERSE DENSITY PROFILE 

We also deduce expression for water density profile along the confinement direction z. 
p{z) is obtained by assuming chemical potential of water fl to be z-dependent and then, a 
partial derivative of ln(Z||) is taken with respect to f3fi{z). At both interfaces i.e., z = 1 
and z = L — 1, the modified fugacity provides additional correction to average density. The 
expression for density profile is given by : 



_ c^(ln(Z||)) . 
d{Pfi{z)) 



Pc is obtained from differentiating Z in Eq. (!28|) . It is the density profile between ideal 
hydrophobic surfaces {us = 0) and is the dominant contribution at all positions. The explicit 
expression for pc{z) upto one- loop order is given in the Appendix ([C]). The interfaces- 
dependent term in Eq. fH3|) can be analyzed using Eq. fl371) . This contribution is only at 
z = 1 and z = L — 1. 

The transverse density profile is shown in Fig.( fT3ll after scaling p{z) with respect to bulk 
density value. At both interfaces there is a characterstic rise in density. From expressions 
of pc{.z) (Appendix [C]) and interface terms (Eq l37|) it is evident that net contribution of 
field correlations is numerically small since density is charge-neutral quantity and linear 0- 
dependent terms tend to cancel each other. Hence, away from interfaces density reaches bulk 
density value rapidly within a distance Many a model simulations in the past computed 
the transverse density profile for water confined between model hydrophobic surfaces. The 



short distance density increase is generically observed 42|, |43| . At ambient conditions the 



magnitude of interfacial density is seen to be typically 1.3 times bulk density value near 
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surfaces with alkane headgroups j43| and independent of L. In our model study we see 
an L- independent increase of magnitude 1.2 for an ideal hydrophobic surface. The under- 
estimation could possibly be due to discrete orientational freedom envisaged in our model. 
Also, alkane head-groups in simulations may have an extra entropy due to fluctuating short 
length polymer chains. 

The rise in interfacial density is also seen for water in the vicinity of hydrophilic surfaces 
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44( 1 . In our model study p{z) between hydrophilic surfaces also displays qualitatively 
similar profile and a lower magnitude of interfacial density compared to that near an ideal 
hydrophobic surface. In all cases the phenomenon is seen to be a consequence of the fact 
that water density has to vanish on the surface. This is compensated by an increase at the 
interface and the system comes back to its bulk equilibrium density within a distance 
from the interface. 

We also calculate density correlations within the interfacial plane and between sites on 
interface and away from interface. Density correlations between any two sites r and r' can 
be calculated from : 

{W{r)W{r')) = ( f-;\^^-v\ ) (44) 

\Zsite{r) Zsite{r)/ 

where, Zsite{f) is site functional at r. To compute density correlations on same interface, 
the site functional at both sites is given by Eq. fl27|) . For density correlations between a site 
on interface and another, away from interface, the site functionals are given by Eqs. (!27l [T6|) 
respectively. //{...} refers to the term proportional to fi in the respective site functional. 
The connected part of the correlation is given by < W{r)W{r') >c = (< W{r)W{r') > — < 
W{r) >< W{r ) >). The explicit expression in each context is deduced upto one-loop order 
in terms of Qri, Q<t> and are given in Appendix ( iDl) . 

Density correlations scaled appropriately with respect to bulk density value are plotted 
in Fig. lfT^ . The plot corresponds to h = 3.58. The figure essentially indicates density 
correlations do not extend beyond few molecular diameters from the interface. Also, there 
is no significant difference between correlations within an interface and that of between 
interface and non-interface sites. 

Similarly, orientational correlations can also be analyzed using the expressions for ori- 
entational weights given in Appendix (|B]). Their effect persists upto longer distance away 
from interface, proportional to the long correlation length of field. 
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Appendix A: Correlation functions in bulk water 



The precise expressions for propagator functions P^^, P^^ deduced in MMF theory are : 



9 



10 



— - /X + A -— + n - -r - 77 



9 



V A' 



10 



+ 4A' + 16/i' - {u' + 2A' - 
^A(l-A) + .' 



(Al) 



(A2) 



1 3 

where, A = - y^(l — cos(/cj)) and v — 2WZo, A' = A/Zq, — 90iJ,/Zo are scaled fugacities. 



They vary between and 1. From zeroth order partition function, DB = u , HB = A , 

The Green's functions for rj, fields in bulk water are given by : 

} dk exp (ik ■ in - 

' ^ (A3) 



and similarly for field. 7^, 7^ are position indices for any two sites. 
For large r — \i\ — 7^|, Qr^if) is of the functional form, 

r ( \ exp {-rjQ . 

where, to the leading order, 

{in)-' 



3 


9 
10 




P 




3 






8 




9 

10 




P 



50/9 



sin I -tan ^\ — I -prz — p 



27 V25 



cos|itan-^/|(|-p 



(A4) 

(A5) 
(A6) 



The form of implies periodic peaks whose amplitudes fall-off exponentially with distance. 
The Qfj, correlator, in addition to oscillatory behavior at short distances, takes the following 
asymptotic form for large r : 



exp 



where. 



-1 



i 



6 



. 5 DB ^ 



6 



1 + ^(4 - /i) - 1 



(A7) 



(A8) 



The above expression for is given to the leading order. 
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Appendix B: Orientational weight C{r],(p) 



The orientational weight for water state in bulk water is given by : 



exp 



He, = 0, ±1 
a = ±1,±2,±3 



Y^{Hl{r)r]{r + e„) + H^{r)(j){r + e„)) 



where, prime indicates summation is subject to constraints Eqs. (fTp!) . 

About the mean field configuration t] = (p = fields are expanded upto quadratic order. 
The C{ri, 0) is then given by : 



C(r/,0) ~ 90 



2^ — V 1 — \ 2 2 — \ 1 — \ 2 2 — \ 

^ + y 2^^" - 3 Z^^- - 5 2^^°^/^ - 3 Z^*^" + 

a a a, 13 a a,f3 



(Bl) 



where rja = ri{r + e^) and = 0(r + Ca) ■ 

The orientational weights for affected orientations of an interfacial water near hydrophobic 
surface are denoted by C'(?7,0). With boundary condition rj = (f) = on surface sites, it is 
given by : 



C'{ri,(j)) ~ 60 



3i ^ — \ 3 — \ 2 3 — ^ 3 — ^ 2 1 \^ — ^ 



a,/3 



(B2) 



Near a hydrophilic surface. 



C'^ivA) - 30 



a,/3 



a, {3 a 



(B3) 



Appendix C: Density profile 



Pc{z) in Eq. (H3!) is given by. 



d{f3j2{z)) 



-y ^i)Ti + {-\'fl')T2 + - /i')T3 



(CI) 
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where, 

Ti = (l-2{u' + 2\'))g,{r,r) + ^J2^,ir,r + e^) + g^ir,r) (C2a) 

a 

T2 = + 2X'))g,{r,r) + ^J2g,{r,r + e^) (C2b) 

a 

4 ' 4 4 

^3 = — (z/' + 2A') y^g^(r, r + e^) + — gr?(r + e», r + Ca) + - ^^(r + e^, r + e^) 

Q; ct a 

+ ^ (1 - ^) 5Z ^^(^ + ^ + ^a') - ^ 5^ + e,, r- + e,. ) (C2c) 

where, r = {x, y, z) is a site position, r + is a near-neighbor site in direction. The 
Green's functions ^,;(ri, r2), ^^(ri, are computed using Eg. (139!) . The scaled fugacities v , 
a', /i' are as defined in Appendix 



Appendix D: Density correlations 



The connected part of density correlation between sites on same interface, to the leading 
order, is given by the expression : 



{W{r)W{r)), ^ -ifi'f 



V + 2ATa,(r, r') + ( ^ fl + ^ - /i') ) E ^^(^ + ^' + 



/i'^^ (i/' + 2A')^6;,(r,r' + e„' 



(Dl) 



where, r, r' are arbitrary sites on same interface; r + Cq,, r' + are respective near- 
neighbor sites in the directions e^, respectively, ^rjl'^i? ''"2)5 ^0('"i5''^2) can be computed 
from Eq.( l39|) . The scaled fugacities z/'. A', /i' are as defined in Appendix (1A|) . 

The density correlation between a site on interface and another, away from interface is 
given by. 



[z.' + 2A')2^,(r,r') + -(l-/x')M , ^ 



2. ' ,S / 3z/5 



v\ 1 



r,r + e„ 



(D2) 
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where, r is any site on interface, r' is away from interface; r + Ca, r' + e^' are their respective 
near-neighbor sites. 

Orientational correlations can also be evaluated using the expression for orientational 
weights given in Appendix ([B]) and can be computed using known expressions for Green's 
functions. 
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FIG. 1. Allowed configurations : a water site with two hydrogen arms (+) and two lone-pair arms 
(-) on links around the site (Eqs ll|2[) . A hydrogen bond occurs when a hydrogen arm (+) and 
a lone-pair arm (-) of two molecules meet at a site, (right bottom corner) unit vectors on cubic 
lattice. 
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FIG. 2. Disallowed configurations : non-zero bond arms of same type of two molecules meeting at 
a site; more than two non-zero arms meeting at a site. 
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FIG. 3. Water confined between macroscopic surfaces. iSi, S2 are surface planes at z = and 
z = L respectively; I\, I2 are their respective interfaces at z = 1 and z = L — 1. 
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FIG. 4. (Color online) Inverse temperature /3 and correlation lengths as a function of h. (3 

is measured in units of hydrogen bond strength A in the model. Correlation lengths are expressed 
in lattice units. 
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FIG. 5. (Color online) Different contributions to Gtot for two hydrophobic surfaces (1/5^ = 1^32 — 
—0.5). The curves are plotted for h = 3.58. In the order from top to bottom the curves correspond 
to 7s (i) — 75(00), Gr, Gc and Gtot{L) — Gtot{oo) respectively. The free energy densities are 
measured per unit hydrogen bond strength. 
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FIG. 6. (Color online) Force between two hydrophobic surfaces {vsi = ^^52 = —0.5). Top (red) 
curve corresponds to h = 3.03, middle (green) curve : h = 3.58, bottom (blue) : h = 3.75. Force 
is measured per unit hydrogen bond strength per unit lattice distance. 
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FIG. 7. (Color online) Gc as a function of L. Top (red) curve corresponds to h = 3.03, middle 
(green) curve : h = 3.58, bottom (blue) : h = 3.75. 
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FIG. 8. (Color online) 75 for a hydrophobic surface (z/g = —0.5). Top (red) curve corresponds to 
h = 3.03, middle (green) curve : h = 3.30, bottom (blue) : h = 3.58. 
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FIG. 9. (Color online) Gr for two hydrophobic surfaces {vsi = 1^82 = —0.5). For L > 6, top (red) 
curve corresponds to h = 3.03, middle (green) curve : h = 3.58, bottom (blue) : h = 3.75. 
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FIG. 10. (Color online) Gr for hydrophilic surfaces {vsi = ^82 — —0.9). (+,+) curve corresponds 
to similar type of hydrophilic surfaces and (+, — ), to dissimilar type. Both curves are plotted for 
h = 3.58. 
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FIG. 11. (Color online) Force between two hydrophilic surfaces (ys^ = — —0.9). (+,+) 
indicates similar type of hydrophilic surfaces and (+, — ) indicates dissimilar type. Both curves 
correspond to h = 3.58. 
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FIG. 12. (Color online) Force between a hydrophilic (z/gj = —0.9) and a hydrophobic surface 
{i'S2 = —0.5). Top (red) curve corresponds to h = 3.03, middle (green) curve : h = 3.58, bottom 
(blue) : h = 3.75. (±, 0) indicates that force is between a positively (negatively) charged hydrophilic 
surface and a hydrophobic surface. 
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FIG. 13. (Color online) Transverse density profile for water between two hydrophobic surfaces 
separated by distance L = 16 and (inset) L = 6. Here, h = 3.58. The steeper (red) curve 
corresponds to ideal case z/^^ = z/^j = and the other (green) curve corresponds to vs^ = 1^82 = 
-0.5. 
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FIG. 14. (Color online) Density correlations near a hydrophobic interface {vs = —0.5) scaled 
appropriately with respect to bulk density value aX h = 3.58. Correlations are between a reference 
site on interface {z = 1) and an arbitrary site on a plane defined by its z coordinate. Distance 
between the two sites is measured using Euclidean metric. 
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